
clear
set more off

scalar Do_counterfactuals=1

*********************************************************
capture program drop prog_inputdata
program define prog_inputdata
	args folder inputcsv var type

	local T=15
	local Tplus1=16
	local Tminus1=14
	
	clear
	insheet using "output-`folder'/sim`inputcsv'.csv"
	
	if (`type'==0) {
		drop v`Tplus1' 
		forval t=1/`T' {
			rename v`t' `var'`t'
			cap replace `var'`t'="0" if `var'`t'=="-Inf" | `var'`t'=="NaN"
			cap destring `var'`t', replace
		}
	}
	else if (`type'==1) {
		egen `var'=rowmax(v1 v2 v3 v4 v5 v6 v7 v8 v9 v10 v11 v12 v13 v14 v15)
		keep `var'
		replace `var'=2-`var'
	}
	else if (`type'==2) {
		keep v1
		rename v1 `var'
		destring `var', replace force

	}
	else if (`type'==3) {
		drop v`Tplus1' 
		forval t=1/`T' {
			rename v`t' `var'`t'
			cap replace `var'`t'="1000000" if `var'`t'=="Inf" | `var'`t'=="NaN"
			cap destring `var'`t', replace
		}
	}
	else if (`type'==4) {
		keep if v1==1 //only keep those starting at age 1
		drop v17
		forval t=2/`Tplus1' {
			local s=`t'-1
			rename v`t' `var'`s'
			cap replace `var'`s'="0" if `var'`s'=="-Inf" | `var'`s'=="NaN"
			cap destring `var'`s', replace
		}			
	}

	gen person=_n

	save "output-`folder'/sim`inputcsv'.dta", replace
end
*********************************************************
qui prog_inputdata base health_fromestimated health 4

qui prog_inputdata base hlth health 0
qui prog_inputdata base marstat married 0
qui prog_inputdata base incp incp 0
qui prog_inputdata base inck inck 0
qui prog_inputdata base inctypek inctypek 2
qui prog_inputdata base socinsp SIp 0
qui prog_inputdata base socinsk SIk 0
qui prog_inputdata base conp conp 0
qui prog_inputdata base conk conk 0
qui prog_inputdata base beq beq 0
qui prog_inputdata base mu mu 0
qui prog_inputdata base labor labor 0
qui prog_inputdata base family family 0
qui prog_inputdata base kap kappa 0
qui prog_inputdata base insd insurance 1
qui prog_inputdata base wstarmar wstar 2
qui prog_inputdata base wstardiv wstardiv 2
qui prog_inputdata base assp assetp 3
qui prog_inputdata base assk assetk 3
qui prog_inputdata base savp savp 3
qui prog_inputdata base savk savk 3

if(Do_counterfactuals==1) {
	qui prog_inputdata cash insd insurance_cash 1
	qui prog_inputdata cash wstarmar wstar_cash 2	
	qui prog_inputdata cash socinsp SIp_cash 0
	qui prog_inputdata cash family family_cash 0
	qui prog_inputdata cash labor labor_cash 0
	qui prog_inputdata cash conp conp_cash 0
	qui prog_inputdata cash conk conk_cash 0	
	
	qui prog_inputdata cash30 insd insurance_cash30 1
	qui prog_inputdata cash30 wstarmar wstar_cash30 2	
	qui prog_inputdata cash30 family family_cash30 0
		
	qui prog_inputdata cash50 insd insurance_cash50 1
	qui prog_inputdata cash50 wstarmar wstar_cash50 2		
	qui prog_inputdata cash50 family family_cash50 0
	
	qui prog_inputdata noinf insd insurance_noinf 1
	qui prog_inputdata noinf socinsp SIp_noinf 0
	qui prog_inputdata noinf labor labor_noinf 0
	qui prog_inputdata noinf conp conp_noinf 0
	qui prog_inputdata noinf conk conk_noinf 0
}

**********************************************************************************************
*TABLE IN PAPER ON HEALTH STATUS PROBABILITIES 
*from estimated probabilities, not true transitions
use "output-base/simhealth_fromestimated.dta", clear

egen everhealthy=anymatch(health1 health2 health3 health4 health5 health6 health7 health8 health9 health10 health11 health12 health13 health14 health15), values(1)
egen everkindasick=anymatch(health1 health2 health3 health4 health5 health6 health7 health8 health9 health10 health11 health12 health13 health14 health15), values(2)
egen eververysick=anymatch(health1 health2 health3 health4 health5 health6 health7 health8 health9 health10 health11 health12 health13 health14 health15), values(3)
egen evereithersick=anymatch(health1 health2 health3 health4 health5 health6 health7 health8 health9 health10 health11 health12 health13 health14 health15), values(2 3)
egen wavesdead=anycount(health1 health2 health3 health4 health5 health6 health7 health8 health9 health10 health11 health12 health13 health14 health15), values(0)
gen lifeexpect_65=(30-2*wavesdead-1)
gen lifeexpect_75=(20-2*wavesdead-1) if health6!=0

gen healthy_65=(health1==1)
gen kindasick_65=(health1==2)
gen verysick_65=(health1==3)
gen eithersick_65=(inlist(health1,2,3))
gen dead_65=(health1==0)

gen healthy_75=(health6==1)
gen kindasick_75=(health6==2)
gen verysick_75=(health6==3)
gen eithersick_75=(inlist(health6,2,3))
gen dead_75=(health6==0)


gen healthy_85=(health11==1)
gen kindasick_85=(health11==2)
gen verysick_85=(health11==3)
gen eithersick_85=(inlist(health11,2,3))
gen dead_85=(health11==0)

gen healthy_93=(health15==1)
gen kindasick_93=(health15==2)
gen verysick_93=(health15==3)
gen eithersick_93=(inlist(health15,2,3))
gen dead_93=(health15==0)

*TABLE:
sum lifeexpect* dead_* everhealthy everkindasick eververysick evereithersick
sum healthy_65 kindasick_65 verysick_65 eithersick_65
sum healthy_75 kindasick_75 verysick_75 eithersick_75 if dead_75!=1
sum healthy_85 kindasick_85 verysick_85 eithersick_85 if dead_85!=1
sum healthy_93 kindasick_93 verysick_93 eithersick_93 if dead_93!=1

*CALCULATING LOAD ON LONG-TERM CARE INSURANCE
reshape long health, i(person) j(age)
gen tempeversick2=(health==2)
bys person: egen eversick2=max(tempeversick2)
bys person: egen numsick2=sum(tempeversick2)
gen tempeversick3=(health==3)
bys person: egen eversick3=max(tempeversick3)
bys person: egen numsick3=sum(tempeversick3)
bys person: gen yrsalive=_n if (inlist(health,1,2,3) & health[_n+1]==0) | _n==15
replace yrsalive=2*yrsalive-1
drop tempeversick*
sum numsick2
scalar avsick2=r(mean)
sum numsick3
scalar avsick3=r(mean)
scalar avbenefit=20000*avsick2+44350*avsick3
disp "average total benefit: " avbenefit
sum yrsalive
scalar avpremium=1/.82*avbenefit/r(mean)
*scalar avpremium=1/1*avbenefit/r(mean)
*scalar avpremium=1/.70*avbenefit/r(mean)
*scalar avpremium=1/.50*avbenefit/r(mean)
disp "average total premiums: " avpremium
disp "average load: " 1-avbenefit/(avpremium*r(mean))

*PREMIUM IF FAIR AND FULL INSURANCE
scalar avbenefit_fairfull=20000*avsick2+61700*avsick3
disp "average fair full total benefit: " avbenefit_fairfull
sum yrsalive
scalar avpremium_fairfull=1/1*avbenefit_fairfull/r(mean)
disp "average total premiums fair full: " avpremium_fairfull
disp "average load fair full: " 1-avbenefit_fairfull/(avpremium_fairfull*r(mean))

**********************************************************************************************

**********************************************************************************************
use "output-base/simhlth.dta", clear
merge 1:1 person using "output-base/siminsd.dta", assert(match) nogen
merge 1:1 person using "output-base/simlabor.dta", assert(match) nogen
merge 1:1 person using "output-base/simfamily.dta", assert(match) nogen
merge 1:1 person using "output-base/simwstarmar.dta", assert(match) nogen
merge 1:1 person using "output-base/simassp.dta", assert(match) nogen
merge 1:1 person using "output-base/simassk.dta", assert(match) nogen
merge 1:1 person using "output-base/simincp.dta", assert(match) nogen
merge 1:1 person using "output-base/siminck.dta", assert(match) nogen
merge 1:1 person using "output-base/siminctypek.dta", assert(match) nogen
merge 1:1 person using "output-base/simsocinsp.dta", assert(match) nogen
merge 1:1 person using "output-base/simsocinsk.dta", assert(match) nogen
merge 1:1 person using "output-base/simconp.dta", assert(match) nogen
merge 1:1 person using "output-base/simconk.dta", assert(match) nogen
merge 1:1 person using "output-base/simmarstat.dta", assert(match) nogen
merge 1:1 person using "output-base/simmu.dta", assert(match) nogen
merge 1:1 person using "output-base/simbeq.dta", assert(match) nogen
merge 1:1 person using "output-base/simkap.dta", assert(match) nogen
merge 1:1 person using "output-base/simsavp.dta", assert(match) nogen
merge 1:1 person using "output-base/simsavk.dta", assert(match) nogen

forval x=1/15 {
	cap destring assetp`x', force replace
	cap destring assetk`x', force replace
	cap destring savp`x', force replace
	cap destring savk`x', force replace
}

if (Do_counterfactuals==1 ) {

	merge 1:1 person using "output-cash/siminsd.dta", assert(match) nogen
	merge 1:1 person using "output-cash/simwstarmar.dta", assert(match) nogen
	merge 1:1 person using "output-cash/simsocinsp.dta", assert(match) nogen
	merge 1:1 person using "output-cash/simfamily.dta", assert(match) nogen
	merge 1:1 person using "output-cash/simlabor.dta", assert(match) nogen
	merge 1:1 person using "output-cash/simconp.dta", assert(match) nogen
	merge 1:1 person using "output-cash/simconk.dta", assert(match) nogen
	
	merge 1:1 person using "output-cash30/siminsd.dta", assert(match) nogen
	merge 1:1 person using "output-cash30/simwstarmar.dta", assert(match) nogen
	merge 1:1 person using "output-cash30/simfamily.dta", assert(match) nogen
	
	merge 1:1 person using "output-cash50/siminsd.dta", assert(match) nogen
	merge 1:1 person using "output-cash50/simwstarmar.dta", assert(match) nogen	
	merge 1:1 person using "output-cash50/simfamily.dta", assert(match) nogen

	merge 1:1 person using "output-noinf/siminsd.dta", assert(match) nogen
	merge 1:1 person using "output-noinf/simsocinsp.dta", assert(match) nogen
	merge 1:1 person using "output-noinf/simlabor.dta", assert(match) nogen
	merge 1:1 person using "output-noinf/simconp.dta", assert(match) nogen
	merge 1:1 person using "output-noinf/simconk.dta", assert(match) nogen

	reshape long health labor family incp inck assetp assetk SIp SIk conp conk mu beq kappa savp savk ///
				 SIp_noinf labor_noinf conp_noinf conk_noinf ///
				 SIp_cash family_cash labor_cash conp_cash conk_cash family_cash30 family_cash50, i(person) j(age)
}
else {
	reshape long health labor family incp inck assetp assetk SIp SIk conp conk mu beq kappa savp savk, i(person) j(age)
}

replace age=64+age*2

bys person: gen parentdied=(health==0 & inrange(health[_n-1],1,3))
gen tempageparentdied=age if parentdied==1
bys person: egen ageparentdied=max(tempageparentdied)
bys person: gen parentdead=(age>=ageparentdied)

replace insurance=. if health==0
cap replace insurance_noinf=. if health==0
cap replace insurance_cash=. if health==0
cap replace insurance_cash30=. if health==0
cap replace insurance_cash50=. if health==0
replace family=. if health==0
cap replace family_cash=. if health==0
replace assetp=. if health==0
replace assetk=. if health==0 & assetk==0
replace savp=. if health==0
replace SIp=. if health==0
cap replace SIp_cash=. if health==0
cap replace SIp_noinf=. if health==0
replace mu=. if health==0

replace conp=. if health==0
replace conk=. if health==0 & parentdead!=1 

if (Do_counterfactuals==1 ) {
	cap replace insurance_noinf=. if health==0
	cap replace insurance_mcdhalf=. if health==0
	cap replace insurance_mcdhalf_noinf=. if health==0
	cap replace insurance_mcd150=. if health==0
	cap replace insurance_mcd150_noinf=. if health==0
	cap replace insurance_mcdcash=. if health==0
	cap replace insurance_allcash=. if health==0
	cap replace family_cash=. if health==0
	cap replace family_cash30=. if health==0
	cap replace family_cash50=. if health==0
	cap replace assetp_noinf=. if health==0
	cap replace SIp_noinf=. if health==0
	cap replace conp_noinf=. if health==0
	cap replace conk_noinf=. if health==0 & parentdead!=1
	cap replace conp_cash=. if health==0
	cap replace conk_cash=. if health==0 & parentdead!=1
}

gen temp66=1 if age==66 & health!=0
bys person: egen cohort1=max(temp66)

gen assetboth=assetp+assetk
xtile assbothq=assetboth if age==66, nq(5)
xtile asskq=assetk if age==66, nq(5)

xtile asspq=assetp if age==66, nq(5) 
xtile asspq4=assetp if age==66, nq(4) 
xtile asspq10=assetp if age==66, nq(10) 

gen 	apq=1 if assetp<1130.476 & age==66
replace apq=2 if assetp<39336.48 & age==66 & apq==.
replace apq=3 if assetp<118997.5 & age==66 & apq==.
replace apq=4 if assetp<296660.7 & age==66 & apq==.
replace apq=5 if assetp>=296660.7 & age==66 & apq==. & assetp!=.
bys person: egen init_apq=max(apq)

bys person: egen init_assetptemp=mean(assetp) if age==66
bys person: egen init_assetp=max(init_assetptemp)

gen 	apq2=1 if assetp<1310.807 
replace apq2=2 if assetp<43256.64 & apq2==.
replace apq2=3 if assetp<124264.5 & apq2==.
replace apq2=4 if assetp<305867.2 & apq2==.
replace apq2=5 if assetp>=305867.2 & apq2==.

gen 	ins_data=0.0269316 if apq==1
replace ins_data=0.040744 if apq==2
replace ins_data=0.0598065 if apq==3
replace ins_data=0.0863656 if apq==4
replace ins_data=0.1699518 if apq==5

gen 	inf_data=0.3927813 if apq==1
replace inf_data=0.4819159 if apq==2
replace inf_data=0.6164384 if apq==3
replace inf_data=0.5985915 if apq==4
replace inf_data=0.5391621 if apq==5

gen 	FTdata=.50555556 if apq2==1
replace FTdata=.54545455 if apq2==2
replace FTdata=.58097802 if apq2==3
replace FTdata=.64487633 if apq2==4
replace FTdata=.62940141 if apq2==5

gen 	PTdata=.0912037 if apq2==1
replace PTdata=.09090909 if apq2==2
replace PTdata=.08120233 if apq2==3
replace PTdata=.08878092 if apq2==4
replace PTdata=.09859155 if apq2==5

gen SIp_data=0.2031141

**********************************************************************************************
*CONVERT PARAMETER ESTIMATES INTO DOLLARS

scalar gammaest=2.1143038768435756
scalar formalcare=-2E-7
scalar guilt=-3E-7
scalar medicaidpref=-4.45076144931274663E-006

*interpret formal care preference in terms of consumption: take median consumption of parent 
sum conp if inlist(health,1,2,3), det 
local medconp=r(p50)
disp `medconp'
disp (`medconp'^(1-gammaest))/(1-gammaest)
*what consumption change would you need to get value of formal care but no formal care:
local altconp= (((`medconp'^(1-gammaest))/(1-gammaest) + formalcare)*(1-gammaest))^(1/(1-gammaest)) 
disp `medconp' -`altconp'

*interpret guilt in terms of consumption: take median consumption of child while parent alive
sum conk if inlist(health,1,2,3), det 
local medconk=r(p50)
disp `medconk' 
*what consumption change would you need to get value of guilt but no guilt:
local altconk= (((`medconk'^(1-gammaest))/(1-gammaest) + guilt)*(1-gammaest))^(1/(1-gammaest)) 
*so change in consumption is 
disp `medconk' -`altconk'

*interpret medicaid preference in terms of consumption: take median consumption of parent
sum conp if inlist(health,1,2,3), det 
local medconp=r(p50)
disp `medconp'
*what consumption change would you need to get value of medicaid but no medicaid (formal care in either case):
local altconp= (((`medconp'^(1-gammaest))/(1-gammaest) + medicaidpref)*(1-gammaest))^(1/(1-gammaest)) 
disp `medconp' -`altconp'
**********************************************************************************************

*************************************************************************************
*CONSUMPTION INSURANCE against income and health shocks - baseline and counterfactuals

preserve
bys person: gen laghealth=health[_n-1]
gen logincp=log(incp)
gen loginck=log(inck)
bys person: gen diff_loginck=loginck-loginck[_n-1]
foreach var in conp conk conp_noinf conk_noinf conp_cash conk_cash {
	keep if inrange(`var',0,1000000)
	gen log`var'=log(`var')
	bys person: gen diff_log`var'=log`var'-log`var'[_n-1]
}
local controls ""
eststo clear
eststo: reg diff_logconp i.health `controls' if laghealth==1 & health!=0, cluster(person)
eststo: reg diff_logconp_noinf i.health `controls' if laghealth==1 & health!=0, cluster(person)
eststo: reg diff_logconp_cash i.health `controls' if laghealth==1 & health!=0, cluster(person)
eststo: reg diff_logconk i.health `controls' if laghealth==1 & health!=0, cluster(person)
eststo: reg diff_logconk_noinf i.health `controls' if laghealth==1 & health!=0, cluster(person)
eststo: reg diff_logconk_cash i.health `controls' if laghealth==1 & health!=0, cluster(person)
esttab, se star(* 0.10 ** 0.05 *** 0.01) b(3)
restore

*************************************************************************************

********************************************************************************

*MODEL FIT: MEDICAID
sum SIp SIp_data

*MODEL FIT: LABOR SUPPLY
tab labor if health==1 //healthy
tab labor if inlist(health,2,3) //sick

*MODEL FIT: BEQUESTS 
preserve
drop if health==0
sort person age
by person: gen lastpd=(_n==_N) & age!=94
replace beq=. if lastpd!=1
sum beq, det
sum assetp if lastpd==1, det
pctile ptile_amount_model=assetp if lastpd==1, nq(1000)
gen ptile=_n
merge 1:1 ptile using "..\bequest_percentiles_data.dta"
replace ptile=ptile/10
replace ptile_amount_data=ptile_amount_data/1000
replace ptile_amount_model=ptile_amount_model/1000
twoway line ptile ptile_amount_data if ptile_amount_data<1000, lcolor(black) lpattern(solid) lwidth(thick) ///
    || line ptile ptile_amount_model if ptile_amount_model<1000, lcolor(black) lpattern(dash) lwidth(thick) ///
    legend(order(1 2) label(1 "Data") label(2 "Model")) ///
	xtitle(Bequest ($1,000s)) ytitle(Percentile of the bequest distribution) ///   
	graphregion(color(white)) xsize(4.5)
graph export "bequest_cdf.pdf", as(pdf) replace	
restore

*HOW DOES MU EVOLVE
cap drop mu_*
gen 	mu_pct = 0.0001 if mu==1
replace mu_pct = 0.1 if mu==2
replace mu_pct = 0.2 if mu==3
replace mu_pct = 0.3 if mu==4
replace mu_pct = 0.4 if mu==5
replace mu_pct = 0.5 if mu==6
replace mu_pct = 0.6 if mu==7
replace mu_pct = 0.7 if mu==8
replace mu_pct = 0.8 if mu==9
replace mu_pct = 0.9 if mu==10
replace mu_pct = 0.999 if mu==11

sort person age
bys person: gen mu_change=  (mu!=mu[_n-1]) if mu!=. & mu[_n-1]!=.
bys person: replace mu_change=  (mu!=7) if mu!=. & (health[_n-1]==0 | age==66)
bys person: gen mu_increase=(mu>mu[_n-1]) if mu!=. & mu[_n-1]!=.
bys person: replace mu_increase =  (mu>7) if mu!=. & (health[_n-1]==0 | age==66)
bys person: gen mu_decrease=(mu<mu[_n-1]) if mu!=. & mu[_n-1]!=.
bys person: replace mu_decrease =  (mu<7) if mu!=. & (health[_n-1]==0 | age==66)
bys person: gen mu_difference=  (mu_pct-mu_pct[_n-1]) if mu!=. & mu[_n-1]!=.
bys person: replace mu_difference =  (mu_pct-.6) if mu!=. & (health[_n-1]==0 | age==66)


bys person: gen inck_decrease=(inck<inck[_n-1]) if inck!=. & inck[_n-1]!=.
bys person: gen health_decrease=(health>health[_n-1]) if inrange(health,1,3) & inrange(health[_n-1],1,3)
bys person: gen health_increase=(health<health[_n-1]) if inrange(health,1,3) & inrange(health[_n-1],1,3)

preserve
collapse (mean) mu_pct_mean=mu_pct (sd) mu_pct_sd=mu_pct (mean) mu_change, by(age)
graph twoway line mu_pct_mean age, lcolor(black) lwidth(thick) || bar mu_change age, ///
			 graphregion(color(white)) ylabel(0(.1)0.5) xlabel(65(5)95) xtitle(Parent age) ///
			 legend(label(1 Mean bargaining weight) label(2 Percent w/ weight change))
graph export "bargaining_weight.pdf", as(pdf) replace		   
restore

********************************************************************************

*how much consumption for those at consumption floor?
gen 	earn=inck if labor==3
replace earn=inck/2 if labor==2
replace earn=0 if labor==1
table SIp, c(median conp median assetp median assetk median assetboth median earn)

*fiscal externality 
gen 	inck_pretax=inck/(1-.2) if labor==3
replace inck_pretax=inck/(1-.2)/2 if labor==2
replace inck_pretax=0 if labor==1
gen 	inck_pretax_noinf=inck/(1-.2) if labor_noinf==3
replace inck_pretax_noinf=inck/(1-.2)/2 if labor_noinf==2		
replace inck_pretax_noinf=0 if labor_noinf==1
gen 	inck_pretax_cash=inck/(1-.2) if labor_cash==3
replace inck_pretax_cash=inck/(1-.2)/2 if labor_cash==2		
replace inck_pretax_cash=0 if labor_cash==1
foreach var in inck_pretax inck_pretax_noinf inck_pretax_cash {
		replace `var'=`var'/2 //because income is biennial -- to make tax stuff annual
}
gen tax_base=0.2*inck_pretax
gen tax_noinf=0.2*inck_pretax_noinf
gen tax_cash=0.2*inck_pretax_cash
gen tax_diff_noinf=tax_noinf-tax_base
gen tax_diff_cash=tax_cash-tax_base
preserve
collapse (mean) tax_base tax_noinf tax_cash tax_diff_noinf tax_diff_cash if health!=0, by(apq2)
twoway line tax_diff_noinf apq2, lcolor(blue) lwidth(thick) || line tax_diff_cash apq2, lcolor(midgreen) lwidth(thick) ///
	   legend(label(1 "No family care avail.") label(2 "Cash")) yline(0, lwidth(medthick) lcolor(black) lpattern(dash)) ///
	    xtitle(Parent wealth quintile) ytitle(Average difference in taxes relative to benchmark) graphregion(color(white)) xsize(4.5)
graph export tax_differences.pdf, as(pdf) replace	
save "tax_differences_byquintile.dta", replace
restore 
*add in medicaid savings for cash version:
preserve
gen cost_medicaid=0
replace cost_medicaid=20000*2      if health==2 & SIp==1 & family==0 & inlist(insurance,0,2)
replace cost_medicaid=61700*2  if health==3 & SIp==1 & family==0 & inlist(insurance,0,2) 
replace cost_medicaid=17350*2  if health==3 & SIp==1 & family==0 & insurance==1 
gen 	cost_medicaid_cash=0
replace cost_medicaid_cash=20000*2      if health==2 & SIp_cash==1  & family_cash==0 & inlist(insurance_cash,0,2)
replace cost_medicaid_cash=61700*2  if health==3 & SIp_cash==1  & family_cash==0 & inlist(insurance_cash,0,2) 
replace cost_medicaid_cash=17350*2  if health==3 & SIp_cash==1  & family_cash==0 & insurance_cash==1 
gen cost_medicaid_diff_cash = cost_medicaid_cash-cost_medicaid
gen diff_tax_mcd = tax_diff_cash - cost_medicaid_diff_cash
tab apq2, sum(tax_diff_cash)
tab apq2, sum(cost_medicaid_diff_cash)
tab apq2, sum(diff_tax_mcd)
collapse (mean) diff_tax_mcd if health!=0, by(apq2)
save "taxandmedicaid_difference_byquintile.dta", replace
restore

*profits
preserve

keep if cohort1==1
bys person: egen yearsparentdead=total(parentdead)
gen yearsalive=30-2*yearsparentdead

cap drop insurer*
gen 	insurer_prem_temp = 3099.375*2 if insurance==1
gen 	insurer_pay_temp = 20000*2 if insurance==1 & health==2 & family==0
replace insurer_pay_temp = 44350*2 if insurance==1 & health==3 & family==0
egen insurer_prem = total(insurer_prem_temp)
egen insurer_pay = total(insurer_pay_temp)
gen insurer_profit = insurer_prem - insurer_pay
count if insurance==1 & age==66
scalar n_insurance=r(N)
foreach var in prem pay profit {
	sum insurer_`var'
	scalar insurer_`var'2=r(mean)
}
sum yearsalive if insurance==1 & age==66
scalar yearsaliveavg=r(mean)
disp insurer_profit2/n_insurance/yearsaliveavg
disp insurer_prem2/n_insurance/yearsaliveavg
disp insurer_pay2/n_insurance/yearsaliveavg
disp yearsaliveavg
gen load_benchmark= 1-insurer_pay/insurer_prem

drop insurer*
gen 	insurer_noinf_prem_temp = 3099.375*2 if insurance_noinf==1
gen 	insurer_noinf_pay_temp = 20000*2 if insurance_noinf==1 & health==2 
replace insurer_noinf_pay_temp = 44350*2 if insurance_noinf==1 & health==3 
egen insurer_noinf_prem = total(insurer_noinf_prem_temp)
egen insurer_noinf_pay = total(insurer_noinf_pay_temp)
gen insurer_noinf_profit = insurer_noinf_prem - insurer_noinf_pay
count if insurance_noinf==1 & age==66
scalar n_insurance_noinf=r(N)
foreach var in prem pay profit {
	sum insurer_noinf_`var'
	scalar insurer_noinf_`var'2=r(mean)
}
sum yearsalive if insurance_noinf==1 & age==66
scalar yearsaliveavg=r(mean)
disp insurer_noinf_profit2/n_insurance_noinf/yearsaliveavg
disp insurer_noinf_prem2/n_insurance_noinf/yearsaliveavg
disp insurer_noinf_pay2/n_insurance_noinf/yearsaliveavg
disp yearsaliveavg
gen load_noinf=  1-insurer_noinf_pay/insurer_noinf_prem


foreach var in cash cash30 cash50 {
	cap drop insurer*
	gen 	insurer_`var'_prem_temp = 3099.375*2 if insurance_`var'==1 & family_`var'==0
	gen 	insurer_`var'_pay_temp = 20000*2 if insurance_`var'==1 & health==2 & family_`var'==0
	replace insurer_`var'_pay_temp = 44350*2 if insurance_`var'==1 & health==3 & family_`var'==0
	egen insurer_`var'_prem = total(insurer_`var'_prem_temp)
	egen insurer_`var'_pay = total(insurer_`var'_pay_temp)
	gen insurer_`var'_profit = insurer_`var'_prem - insurer_`var'_pay
	count if insurance_`var'==1 & age==66
	scalar n_insurance_`var'=r(N)
	foreach var2 in prem pay profit {
		sum insurer_`var'_`var2'
		scalar insurer_`var'_`var2'2=r(mean)
	}
	sum yearsalive if insurance_`var'==1 & age==66
	scalar yearsaliveavg=r(mean)
	disp insurer_`var'_profit2/n_insurance_`var'/yearsaliveavg
	disp insurer_`var'_prem2/n_insurance_`var'/yearsaliveavg
	disp insurer_`var'_pay2/n_insurance_`var'/yearsaliveavg
	disp yearsaliveavg
	gen load_`var'=  1-insurer_`var'_pay/insurer_`var'_prem	
}

keep load_benchmark load_noinf load_cash load_cash30 load_cash50
sum load*
restore


*inter-vivos transfers
gen 	cohk=(assetk+inck) if labor==3
replace cohk=(assetk+0.5*inck) if labor==2
replace cohk=(assetk+0.0*inck) if labor==1
replace cohk = 23381 if cohk<23381
gen transfer_ktop = (cohk-conk-savk)/2
gen transfer_ptok=-transfer_ktop

gen sick=inlist(health,2,3)
by person: egen bequest=max(beq)
by person: replace bequest=. if _n!=1

eststo clear
eststo: reg transfer_ptok incp assetp inck assetk sick family insurance age if transfer_ktop<500000 & transfer_ktop>-500000 & assetp<10000000 & assetk<10000000 & health!=0, cluster(person)
estadd ysumm
esttab using "transfers_regression.tex", replace booktabs nonotes b(a2) se(a2) ///
			 star(* 0.10 ** 0.05 *** 0.01) stats(N ymean, fmt(a3 a3) label("N" "Mean net transfer")) mtitle("Net transfer (parent to child)") ///
			 varlabel(age "Age of parent" insurance "LTC insurance" inck "Child income" incp "Parent income" ///
			 assetp "Parent assets" assetk "Child assets" sick "Parent needs LTC" family "Informal care" _cons "Constant") 
		
*threat points
*fraction of changes that are an increase:
gen mu_increase_fraction=mu_increase if mu_change==1
gen mu_decrease_fraction=mu_decrease if mu_change==1
sum mu_increase_fraction

eststo clear
eststo: reg mu_increase health_decrease health_increase inck_decrease if inrange(health,1,3), cluster(person)
eststo: reg mu_increase_fraction health_decrease health_increase inck_decrease if inrange(health,1,3), cluster(person)
eststo: reg mu_decrease health_decrease health_increase inck_decrease if inrange(health,1,3), cluster(person)
eststo: reg mu_decrease_fraction health_decrease health_increase inck_decrease if inrange(health,1,3), cluster(person)
esttab using "mu_change_regression.tex", replace booktabs nonotes b(3) se(3) star(* 0.10 ** 0.05 *** 0.01) stats(N, fmt(a3)) ///
	mtitles("Parent threat" "Parent threat, cond." "Child threat" "Child threat, cond.") ///
	varlabel(health_decrease "Parent health deteriorated" health_increase "Parent health improved" inck_decrease "Child income decreased" _cons "Constant")

*what happens when parent gets sick
*find first time parent gets sick, label that event time zero
cap drop agesick firstsick eventtime maxeventtime mineventtime
gen agesick=age if inlist(health,2,3)
bys person: egen firstsick=min(agesick)
gen eventtime=age-firstsick if health!=0
bys person: egen maxeventtime=max(eventtime)
bys person: egen mineventtime=min(eventtime)
gen FT=(labor==3)
gen PT=(labor==2)
gen working=(inlist(labor,2,3))

preserve
local t=4
keep if inrange(eventtime,-`t',`t')
gen highassets=1 if inlist(init_apq,4,5)
replace highassets=0 if inlist(init_apq,1,2)
collapse (mean) FT PT working mu_decrease mu_increase family (median) conp assetp, by(eventtime highassets)
label define highassetsval 0 "Low parent assets" 1 "High parent assets"
label values highassets highassetsval
tab highassets
twoway connected working family mu_decrease mu_increase eventtime, by(highassets, note("") graphregion(color(white))) ///
	color(black black orange orange) lpattern(solid dash solid dash) ///
	lwidth(medthick medthick medthick medthick) msize(small small small small) ///
	legend(label(1 "Working") label(2 "Informal Care") label(3 "Pareto weight decrease") label(4 "Pareto weight increase")) ///
	xline(0) xtitle(Years since parent health shock) graphregion(color(white))
graph export healthshock_outcomes_byassets.pdf, as(pdf) replace
restore

preserve
local t=4
keep if inrange(eventtime,-`t',`t')
collapse (mean) FT PT working mu_decrease mu_increase family (median) conp assetp, by(eventtime insurance)
label define insuranceval 0 "No insurance" 1 "Insurance"
label values insurance insuranceval
twoway connected working family mu_decrease mu_increase eventtime, by(insurance, note("") graphregion(color(white))) ///
	color(black black orange orange) lpattern(solid dash solid dash) ///
	lwidth(medthick medthick medthick medthick) msize(small small small small) ///
	legend(label(1 "Working") label(2 "Informal Care") label(3 "Pareto weight decrease") label(4 "Pareto weight increase")) ///
	xline(0) xtitle(Years since parent health shock) graphregion(color(white))
graph export healthshock_outcomes_byinsurance.pdf, as(pdf) replace
restore

*full time by wealth rather than cohort
gen workingdata=FTdata+PTdata
preserve
label define assetlabel 1 "Wealth Q1" 2 "Wealth Q2" 3 "Wealth Q3" 4 "Wealth Q4" 5 "Wealth Q5"
label values apq2 assetlabel
collapse working workingdata, by(apq2)
graph bar workingdata working, over(apq2) ylabel(0(.2)1.0) yscale(range(0(.2)1)) graphregion(color(white)) xsize(4.5) ///
		  bar(1, color(black)) bar(2, color(gray)) ///
		  ytitle(Full-time work) legend(label(1 "Data") label(2 "Model"))
graph export kidwork_wealthquintile.pdf, as(pdf) replace	
restore

*adverse selection
*construct average total lifetime ltc needs by permanent income tercile
gen health2=(health==2)
gen health3=(health==3)
bys person: egen tothealth2=total(health2)
bys person: egen tothealth3=total(health3)
gen tothealthcost = tothealth2*20000+tothealth3*61700
table incp if age==66 & health!=0, c(mean tothealthcost mean insurance)
table incp if age==66 & health==1, c(mean tothealthcost mean insurance)
table inctypek if age==66 & health!=0, c(mean tothealthcost mean insurance mean family)
table inctypek if age==66 & health==1, c(mean tothealthcost mean insurance mean family)

*"nursing home" rates by wealth
gen formalcare=1-family
tab init_apq if inlist(health,3), sum(formalcare)
preserve
keep if inlist(health,3)
collapse formalcare, by(init_apq)
mat formal_data=[.65351812, .60546875, .45675676, .51736111, .5528169]' //from data preparation files
svmat formal_data
rename formal_data1 formal_data
replace formal_data=. if _n>5
gen data_apq=_n if _n<=5
twoway line formal_data data_apq, color(black) lwidth(thick) lpattern(solid) || line formalcare init_apq, color(black) lwidth(thick) lpattern(dash) ///
	   xtitle(Parent wealth quintile) ytitle("Formal care rate, parents with heavy LTC needs") legend(row(1) label(1 Data) label(2 Model)) ylabel(0(.2)1.0) yscale(range(0(.2)1) noextend) graphregion(color(white)) xsize(4.5)
graph export formalcare_quintile_health3.pdf, as(pdf) replace		   
restore

*******************************************************************************
*WILLINGNESS TO PAY*
/*
wstar	  		: (interpolated asset gridpoint) of money you would have to give 
				  parent to be indifferent between insurance and no insurance
wstar_cash		: same idea, between cash insurance and no insurance

WTP and WTP_cash: convert wstar and wstar_cash to money value

welfare_cash	: difference between WTP_cash and WTP
*/

*convert wstar gridpoint to money value (asset nodes are spaced assetspace apart)
scalar assetspace=1000000/19 //if use 20 gridpoints
cap drop WTP* 
cap drop welfare*
cap drop wstar_assetgrid*

gen wstar_assetgrid=assetspace*(wstar-1)
gen wstar_assetgrid_cash=assetspace*(wstar_cash-1)
gen wstar_assetgrid_cash30=assetspace*(wstar_cash30-1)
gen wstar_assetgrid_cash50=assetspace*(wstar_cash50-1)

gen 	welfare_ik=wstar_assetgrid-assetp if insurance==1
replace welfare_ik=0 if insurance==0

gen 	welfare_cash=wstar_assetgrid_cash-assetp if insurance_cash==1
replace welfare_cash=0 if insurance_cash==0

gen 	welfare_cash30=wstar_assetgrid_cash30-assetp if insurance_cash30==1
replace welfare_cash30=0 if insurance_cash30==0

gen 	welfare_cash50=wstar_assetgrid_cash50-assetp if insurance_cash50==1
replace welfare_cash50=0 if insurance_cash50==0


gen 	welfare_diff=max(welfare_cash,0)-max(welfare_ik,0)
gen 	welfare_diff30=max(welfare_cash30,0)-max(welfare_ik,0)
gen 	welfare_diff50=max(welfare_cash50,0)-max(welfare_ik,0)

foreach var in diff diff30 diff50 {
	replace welfare_`var'=0 if welfare_`var'<0 & insurance!=1
	replace welfare_`var'=. if welfare_`var'<-100000 & insurance==1
}

preserve
keep if cohort1==1 & age==66 //remember can only do insurance stuff with COHORT 1
collapse welfare_diff welfare_diff30 welfare_diff50, by(apq)
twoway connected welfare_diff apq, color(midgreen) lwidth(thick) yline(0, lwidth(medthick) lcolor(black) lpattern(dash)) || ///
	   connected welfare_diff30 apq , color(midgreen) lwidth(thick) lpattern(dash) ||  ///
	   connected welfare_diff50 apq , color(midgreen) lwidth(thick) lpattern(shortdash)  ///
	   graphregion(color(white)) xtitle(Parent wealth quintile) ytitle(Family welfare gain ($)) ///
	   legend(label(1 18% load) label(2 30% load) label(3 50% load) rows(1)) xsize(4.75)   
graph export "welfare_cash.pdf", as(pdf) replace

cap drop apq2
gen apq2=apq
merge 1:1 apq2 using "tax_differences_byquintile.dta", nogen
merge 1:1 apq2 using "taxandmedicaid_difference_byquintile.dta", nogen
gen mvpf=-welfare_diff/tax_diff_cash
gen mvpf_taxandmedicaid=-welfare_diff/diff_tax_mcd 
gen mvpf_taxandmedicaid_neg150 = mvpf_taxandmedicaid
replace mvpf_taxandmedicaid_neg150=160 if mvpf_taxandmedicaid<0

twoway connected mvpf apq, color(midgreen) lwidth(thick) ///
	|| connected mvpf_taxandmedicaid_neg150 apq if apq>=3, color(dkgreen) lwidth(medthick) lpattern(dash) ///
		text(167 1 "({&infinity})", color(dkgreen*.3)) text(167 2 "({&infinity})", color(dkgreen*.3)) ///
		text(167 3 "({&infinity})", color(dkgreen*.3)) ///
	|| connected mvpf_taxandmedicaid_neg150 apq if apq<=3, color(dkgreen*.3) lwidth(medthick) lpattern(dash) ///	
	   graphregion(color(white)) xtitle(Parent wealth quintile) ytitle("MVPF, cash counterfactual") ///
	   legend(order(1 "Taxes only" 2 "Taxes and Medicaid") rows(1)) xsize(4.75)  
graph export "mvpf_cash_all.pdf", as(pdf) replace
	  
restore

********************************************************************************
	   
*MEDICAID USE
preserve

keep if cohort1==1 & health!=0

collapse (mean) SIp SIp_noinf SIp_cash, by(init_apq)
mat mcd_data=[.5571302, .26480836, .14198871, .04321521, .01806452]'
svmat mcd_data
twoway line mcd_data init_apq, lcolor(black) lwidth(thick) || line SIp init_apq, lcolor(black) lwidth(thick) lpattern(dash)  ///
	   legend(order(1 2) label(1 "Data") label(2 "Model")) ///	   
	   xtitle(Parent wealth quintile) ytitle(Medicaid rate) ylabel(0(.2)1) yscale(range(0(.2)1) noextend) ///
	   graphregion(color(white)) xsize(4.5)
graph export "fit_medicaid_quintile.pdf", as(pdf) replace		 

label var SIp "Benchmark"
label var SIp_noinf "No family care avail."
label var SIp_cash "Cash benefit"
cap graph drop medicaid1
twoway connected SIp init_apq, lwidth(thick) color(red) ///
	|| connected SIp_noinf init_apq, lwidth(thick) color(blue) lpattern(dash) ///
	|| connected SIp_cash init_apq, lwidth(thick) color(midgreen) lpattern(shortdash) ///
		xlabel(1(1)5 1 "Poorest") xtitle(Parent wealth quintile) ytitle(Percent of parents on Medicaid) ///
		graphregion(color(white)) legend(rows(1) rowgap(.5) keygap(.5) symxsize(10)) name(medicaid1) xsize(4.75)
graph export "counterfactual1_medicaid.pdf", as(pdf) replace

restore

*COSTS

keep if cohort1==1 & health!=0

gen cost_insurance=0
replace cost_insurance=20000*2     if health==2          & family==0 & insurance==1 //doesn't matter if on Medicaid cuz Medicaid is secondary payer
replace cost_insurance=44350*2 if health==3          & family==0 & insurance==1 //doesn't matter if on Medicaid cuz Medicaid is secondary payer
gen cost_OOP=0
replace cost_OOP=20000*2           if health==2 & SIp==0 & family==0 & inlist(insurance,0,2)
replace cost_OOP=61700*2       if health==3 & SIp==0 & family==0 & inlist(insurance,0,2)
replace cost_OOP=17350*2       if health==3 & SIp==0 & family==0 & insurance==1
gen cost_medicaid=0
replace cost_medicaid=20000*2      if health==2 & SIp==1 & family==0 & inlist(insurance,0,2)
replace cost_medicaid=61700*2  if health==3 & SIp==1 & family==0 & inlist(insurance,0,2) //include 7800 in 'consumption'?
replace cost_medicaid=17350*2  if health==3 & SIp==1 & family==0 & insurance==1 //include 7800 in 'consumption'?
gen cost_family=0
replace cost_family=inck/2*2	   if health==2 		 & family==1
replace cost_family=inck*2	   if health==3 		 & family==1
*add 'ssi'-type part of medicaid? (ie if healthy?)
foreach var in medicaid OOP family insurance {
	sum cost_`var'
	scalar totalcost_`var'=r(N)*r(mean)
	disp totalcost_`var'
}
scalar totalcost = totalcost_insurance + totalcost_OOP + totalcost_family + totalcost_medicaid
disp totalcost
foreach var in medicaid OOP family insurance {
	disp totalcost_`var'/totalcost
}


*NO INFORMAL COUNTERFACTUAL
gen		cost_insurance_noinf=0
replace cost_insurance_noinf=20000*2 	 if health==2 & insurance_noinf==1 //doesn't matter if on Medicaid cuz Medicaid is secondary payer
replace cost_insurance_noinf=44350*2 if health==3 & insurance_noinf==1 //doesn't matter if on Medicaid cuz Medicaid is secondary payer
gen 	cost_OOP_noinf=0
replace cost_OOP_noinf=20000*2           if health==2 & SIp_noinf==0 & inlist(insurance_noinf,0,2)
replace cost_OOP_noinf=61700*2       if health==3 & SIp_noinf==0 & inlist(insurance_noinf,0,2)
replace cost_OOP_noinf=17350*2       if health==3 & SIp_noinf==0 & insurance_noinf==1
gen 	cost_medicaid_noinf=0
replace cost_medicaid_noinf=20000*2      if health==2 & SIp_noinf==1 & inlist(insurance_noinf,0,2)
replace cost_medicaid_noinf=61700*2  if health==3 & SIp_noinf==1 & inlist(insurance_noinf,0,2) //include 7800 in 'consumption'?
replace cost_medicaid_noinf=17350*2  if health==3 & SIp_noinf==1 & insurance_noinf==1 //include 7800 in 'consumption'?
foreach var in medicaid OOP insurance {
	sum cost_`var'_noinf
	scalar totalcost_`var'_noinf=r(N)*r(mean)
	disp totalcost_`var'_noinf
}
scalar totalcost_noinf = totalcost_insurance_noinf + totalcost_OOP_noinf + totalcost_medicaid_noinf
disp totalcost_noinf
foreach var in medicaid OOP insurance {
	disp totalcost_`var'_noinf/totalcost_noinf
}
disp totalcost_medicaid_noinf/totalcost_medicaid
disp totalcost_OOP_noinf/totalcost_OOP
disp totalcost_insurance_noinf/totalcost_insurance
disp totalcost_noinf/totalcost

*CASH COUNTERFACTUAL
gen 	cost_insurance_cash=0
replace cost_insurance_cash=20000*2 	if health==2 & 				  				   insurance_cash==1 //doesn't matter if on Medicaid cuz Medicaid is secondary payer
replace cost_insurance_cash=44350*2 if health==3 & 				  				   insurance_cash==1 //doesn't matter if on Medicaid cuz Medicaid is secondary payer
gen 	cost_OOP_cash=0
replace cost_OOP_cash=20000*2           if health==2 & SIp_cash==0  & family_cash==0 & inlist(insurance_cash,0,2)
replace cost_OOP_cash=61700*2       if health==3 & SIp_cash==0  & family_cash==0 & inlist(insurance_cash,0,2)
replace cost_OOP_cash=17350*2       if health==3 & SIp_cash==0  & family_cash==0 & insurance_cash==1
gen 	cost_family_cash=0
replace cost_family_cash=inck/2*2		if health==2 & 				  family_cash==1 & inlist(insurance_cash,0,2)
replace cost_family_cash=inck*2		if health==3 & 				  family_cash==1 & inlist(insurance_cash,0,2)
replace cost_family_cash=2*max(inck/2-20000,0) if health==2 	& family_cash==1 & insurance_cash==1
replace cost_family_cash=2*max(inck  -44350,0) if health==3 	& family_cash==1 & insurance_cash==1
gen 	cost_medicaid_cash=0
replace cost_medicaid_cash=20000*2      if health==2 & SIp_cash==1  & family_cash==0 & inlist(insurance_cash,0,2)
replace cost_medicaid_cash=61700*2  if health==3 & SIp_cash==1  & family_cash==0 & inlist(insurance_cash,0,2) //include 7800 in 'consumption'?
replace cost_medicaid_cash=17350*2  if health==3 & SIp_cash==1  & family_cash==0 & insurance_cash==1 //include 7800 in 'consumption'?
*add 'ssi'-type part of medicaid? (ie if healthy?)
foreach var in medicaid OOP family insurance {
	sum cost_`var'_cash
	scalar totalcost_`var'_cash=r(N)*r(mean)
	disp totalcost_`var'_cash
}
scalar totalcost_cash = totalcost_insurance_cash + totalcost_OOP_cash + totalcost_family_cash+ totalcost_medicaid_cash
disp totalcost_cash
foreach var in medicaid OOP family insurance {
	disp totalcost_`var'_cash/totalcost_cash
}
disp totalcost_medicaid_cash/totalcost_medicaid
disp totalcost_OOP_cash/totalcost_OOP
disp totalcost_family_cash/totalcost_family
disp totalcost_insurance_cash/totalcost_insurance
disp totalcost_cash/totalcost

bys person: egen personcost_medicaid=total(cost_medicaid)
bys person: egen personcost_medicaid_noinf=total(cost_medicaid_noinf)
bys person: egen personcost_medicaid_cash=total(cost_medicaid_cash)
bys person: replace personcost_medicaid=. if _n!=1
bys person: replace personcost_medicaid_noinf=. if _n!=1
bys person: replace personcost_medicaid_cash=. if _n!=1


mat define pcts=J(2,3,100)
mat pcts[1,2]=100*totalcost_medicaid_noinf/totalcost_medicaid
mat pcts[1,3]=100*totalcost_medicaid_cash/totalcost_medicaid
mat pcts[2,2]=100*totalcost_insurance_noinf/totalcost_insurance
mat pcts[2,3]=100*totalcost_insurance_cash/totalcost_insurance

preserve
clear
svmat pcts
rename pcts1 baseline
rename pcts2 noinf
rename pcts3 cash
gen payer="Medicaid" if _n==1
replace payer="Private insurance" if _n==2
list
graph bar baseline noinf cash if payer=="Medicaid", ///
	  bar(1, color(red) lwidth(thick)) bar(2, lcolor(blue) lwidth(thick) fcolor(blue) fintensity(inten70)) bar(3, lcolor(midgreen) lwidth(thick) fcolor(lime) fintensity(inten30)) outergap(*2) bargap(20) ///
	  yline(100, lpattern(dash) lcolor(black)) ///
	  legend(rows(1) label(1 Benchmark) label(2 No family care avail.) label(3 Cash)) ///
	  ytitle(Total Medicaid spending as % of benchmark) graphregion(color(white)) xsize(4.75) 
graph export "counterfactual_medicaidspending.pdf", as(pdf) replace
restore
